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Simulations of frustrated quantum antiferromagnets suffer from a severe sign problem. We solve 
the ergodicity problem of the loop-cluster algorithm in a natural way and apply a powerful strategy 
to address the sign problem. For the spin | Heisenberg antiferromagnet on a kagome and on a 
frustrated square lattice, a nested cluster algorithm eliminates the sign problem for large systems. 
The method is applicable to general lattice geometries but limited to moderate temperatures. 
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Monte Carlo calculations are a powerful tool for first 
principles investigations of strongly coupled quantum 
systems. Early simulations of quantum spin systems 
were performed with local Metropolis-type algorithms 
They suffered from critical slowing down and were 
thus limited to rather high temperatures. Cluster al- 
gorithms perform non-local updates and are capable of 
substantially reducing critical slowing down. Such al- 
gorithms were first developed by Swendsen and Wang 
for discrete classical Ising and Potts spins Q and then 
generalized by Wolff Q to classical spins with a con- 
tinuous 0{N) symmetry. Improved estimators which 
average over an exponentially large number of configu- 
rations at polynomial cost are an additional benefit of 
cluster algorithms. The first cluster algorithm for the 
spin j quantum Heisenberg model was developed in 
While that algorithm works efficiently only for quantum 
spin chains, the loop-cluster algorithm [5[ is efficient also 
in higher dimensions, and was first applied to the 2-d 
spin J Heisenberg antiferromagnet on a square lattice in 
[6| . The continuous-time variant of the algorithm elim- 
inates the Suzuki- Trotter time-discretization error and 
can reach very low temperatures Q. This method has 
also been used to simulate systems on very large lattices 
[H and with very long correlation lengths [§] . An elegant 
and powerful related method based on stochastic series 
expansion is available as well [lfj . 

Unfortunately, in many cases of physical interest, 
including frustrated quantum spin systems, quantum 
Monte Carlo calculations suffer from a very severe sign 
problem. Using an improved estimator, the sign problem 
of the 2-d classical 0(3) model at vacuum angle 9 = tt has 
been addressed with a variant of the Wolff cluster algo- 
rithm 11]. In that case, some clusters are half-instantons 
also known as merons. Flipping a meron-cluster leads 
to a sign-change of the Boltzmann weight and hence to 
an exact cancellation between two configurations. As a 
consequence, only configurations without meron-clusters 
contribute to the partition function. Restricting the sim- 
ulation to those configurations eliminates the sign prob- 
lem, since all configurations in the zero-meron sector have 
a positive sign. The meron concept has been gener- 
alized to fermionic systems 12| and the meron-cluster 



algorithm has been used to solve a number of very se- 
vere fermion sign problems TH, 14 , [l5| . Unfortunately, 
the meron-cluster algorithm is not generally applicable. 
In fact, as shown in some sign problems are NP- 
complete. Hence, a hypothetical method that can solve 
any sign problem would solve all NP-complete problems 
in polynomial time. This would imply the equality of the 
complexity classes NP=P. Since it is generally believed 
that NP^P, it is expected that a universally applicable 
method that solves all sign problems cannot exist. In this 
paper we construct a nested cluster algorithm which, for 
the first time, is capable of eliminating severe sign prob- 
lems of frustrated antiferromagnets at least at moderate 
temperatures. 

Let us consider the antiferromagnetic spin ^ quantum 
Heisenberg model with the Hamiltonian 
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Here S x is a quantum spin operator located at the site x 
of a lattice A, and J xy > is the antiferromagnetic ex- 
change coupling between a pair of spins located at the 
sites x and y. Although our method can be applied 
directly in the Euclidean time continuum, in order to 
ease its implementation we work in discrete time. De- 
pending on the lattice geometry, the Hamiltonian H = 
Hi + H% + ... + Hm is expressed as a sum of M terms 
Hi which leads to a Suzuki- Trotter decomposition of the 
partition function 

Z=Ttexp(-j3H) 
= ]hnTr[exp(-eH 1 )exp(-eH 2 )...exp(-eH M )] N .(2) 

Here the inverse temperature j3 = 1/T = Ne represents 
the extent of a periodic Euclidean time interval, which 
is divided into N discrete time steps of size e. Each 
Hi is a sum of mutually commuting pair interactions 
h xy = J X yS x • S y on a set of disconnected bonds. Insert- 
ing complete sets of spin states s x ,t — ±5 = ^> i between 
the factors exp(— sHi) in eq.©, the partition function is 
expressed as a path integral over all spin configurations 
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[s] 

Z = ]TSign[ S ]exp(-S[ S ]). (3) 
M 

The weight of a configuration is a product of con- 
tributions from individual space-time plaquettes cor- 
responding to the two-spin transfer matrix elements 
{sx,tSy,t\exp(-eh xy )\s X} t+iSy t t+i)- In the basis | IT), 
tl)> I lt)i I II): U P to an irrelevant overall factor, the 
two-spin transfer matrix takes the form 
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with A = exp(—eJ xy /2) and B = s'mh(eJ xy /2). The 
off-diagonal transfer matrix elements are negative. The 
product of the negative signs over all space-time plaque- 
ttes defines the total Sign[s] = ±1 of a spin configuration. 
The remaining factor exp(— S[s]) represents a positive 
Boltzmann weight which can be interpreted as a proba- 
bility and thus can be used for importance-sampling in a 
Monte Carlo simulation. 

When one samples the system using the positive weight 
cxp(— S[s\), one must include Sign[s] in the measured 
observables 0[s] and expectation values are given by 

(O) = ^E°W Sign[ S ]exp(-S[.s]) = { ° ' Sl g"> + . (5) 
Z ^ (Sign) + 

Here the index + refers to expectation values in the simu- 
lated ensemble with positive Boltzmann weights and par- 
tition function Z + = X)[ s ] ex P(~ S[ S Y) such that 

(Sign)+ = — VSign[s]exp(-5[s]) = — 

z+ w z+ 

~ exp(-A//3F). (6) 

Here V is the spatial volume and A/ is the difference 
between the free energy densities of the original ensem- 
ble with the weight Sign[s] exp(— S[s]) and the simulated 
ensemble with the positive weight exp(— S[s]). The ex- 
pectation value of the sign is exponentially small in the 
space-time volume j3V . Since it is obtained as a Monte 
Carlo average of contributions Sign[s] = ±1, one needs an 
exponentially large statistics in order to accurately mea- 
sure (Sign)-).. This is impossible in practice and gives rise 
to a very severe sign problem. 

How can one increase the statistics by an exponential 
factor without investing more than a polynomial numeri- 
cal effort? The meron-cluster algorithm [ll|, [l2| achieves 
this by constructing an improved estimator for the sign. 
Like the meron-cluster algorithm, the method presented 
here is based on the loop-cluster algorithm Q which deco- 
rates a spin configuration with bonds connecting spins to 



closed loop-clusters. The four spins on a space-time pla- 
quette are connected in pairs. In fact, A and B in eq.((4]) 
represent weights of two possible bond configurations on 
a space-time plaquette. The weight A corresponds to 
bonds connecting the spins s x j and s y j with their time- 
like neighbors s Xj t+i and s y j+i, while B corresponds to 
space-like bonds connecting s Xj t with s v j and s x ^+i with 
s Vl t+i- Sites connected by bonds form a closed oriented 
loop-cluster. Up to an overall spin-flip of the entire clus- 
ter, the spin configuration on a cluster is determined by 
the cluster geometry. Time-like bonds connect parallel 
spins, while space-like bonds connect anti-parallel spins. 
Integrating out the spins, the partition function can be 
expressed as a sum over bond configurations [b] 

Z = ^2Sign[b]A nA B nB 2 Nc , (7) 

lb] 

Here ua is the number of time-like and ns is the number 
of space-like plaquette break-ups, while Nq is the num- 
ber of loop-clusters. The factor 2 Nc arises because each 
cluster has two possible spin orientations. The partition 
function can be sampled by a Metropolis update of the 
plaquette break-ups. Remarkably, while the original clus- 
ter algorithm which operates on spins and bonds never 
changes the sign and is thus not ergodic [l7|, the algo- 
rithm which operates only on bonds (after the spins have 
been integrated out) is ergodic and still avoids unnatural 
freezing. Interestingly, Sign[s] remains invariant under 
cluster flips, i.e. all clusters are non-merons. However, in 
this case the meron-cluster algorithm does not solve the 
sign problem because almost half of the configurations 
in the zero-meron sector have a negative sign [17] . Since 
it does not change under spin flips, Sign[s] = Sign[6] is 
uniquely determined by the bond configuration. It is 
important to note that the sign can be expressed as a 
product of cluster signs Sign[fo] = J| c Sign c . Depend- 
ing on the orientation of a cluster, each space-like break- 
up contributes a factor ±i to the two clusters traversing 
the corresponding space-time plaquette. By construc- 
tion, each cluster traverses an even number of space-like 
break-ups, and hence Sign c = ±1. 

We distinguish space-time plaquettes shared by two 
different clusters from internal plaquettes belonging en- 
tirely to one cluster. Updating the break-up on a space- 
time plaquette shared by two different clusters does not 
lead to a sign-change. Only updates of cluster-internal 
plaquettes may change the sign. We apply the follow- 
ing method to construct an improved estimator for the 
sign. Once a statistically independent bond configura- 
tion has been produced by the cluster algorithm, we per- 
form an inner Monte Carlo simulation by updating only 
the cluster-internal plaquette break-ups. Each cluster C 
defines the set of lattice sites Ac contained in C. The 
inner Monte Carlo algorithm generates clusters with dif- 
ferent orientations that visit all sites of Ac in different 
orders, thus contributing different values of Sign c . In 
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this process, break-ups that lead to the decomposition 
of Ac into separate clusters must be rejected. The inner 
Monte Carlo algorithm estimates an average (Sign c )j for 
each set of sites Ac . Since the different sets are indepen- 
dent, the improved estimator of the sign is given by 

(Sign) 4 =n(Sign c ) i :. (8) 

Ac 

Remarkably, the nesting of an outer and an inner cluster 
algorithm achieves exponential error reduction at poly- 
nomial cost. A similar strategy was very successfully 
applied to the measurement of exponentially suppressed 
Wilson loops in lattice gauge theory [l8[ as well as to 
quantum impurity models [1 91 ] . Correlation functions and 
susceptibilities can also be measured with improved es- 
timators. Let us consider the staggered magnetization 
operator M s = Y^ x z x S x . Here z x is a stagger factor de- 
pending on the sub-lattice to which the site x belongs. 
The corresponding staggered susceptibility 

(A/ s 2 Sign) + ((M s 2 Sign). t ) + 
Xs /3V(Sign) + /3y((Sign) l >+ ' 1 ' 

is obtained from an improved estimator which is given in 
terms of M s = J2c M ^c with M sC = J2(x,t)&c z x$x,t as 

(M 2 Sign) 4 = ^(M^Signc), JJ < Si S n c>- ( 10 ) 

A c A C /^A C 

In which cases will the nested cluster algorithm elim- 
inate or at least substantially reduce the sign problem? 
Since some sign problems are NP-hard, it is expected that 
any method will fail at least in those cases. The nested 
cluster algorithm fails to solve the sign problem when a 
cluster fills almost the entire volume, because then the 
inner Monte Carlo algorithm becomes inefficient. Since 
large clusters necessarily arise in the presence of large 
correlation lengths, the nested cluster algorithm does not 
work efficiently in low-temperature ordered phases. 

Even in the absence of long-range order, cluster al- 
gorithms may become inefficient if the clusters grow to 
unphysically large sizes beyond the physical correlation 
length. This potential problem is prevented when there is 
a reference configuration that limits cluster growth [l5| . 
For the antiferromagnet on the square lattice the refer- 
ence configuration is given by the classical Neel state, i.e. 
all spins in a loop-cluster are in a staggered pattern. The 
cluster-size squared is then tied to the staggered suscep- 
tibility which protects the clusters from growing to un- 
physically large sizes. Also for frustrated systems it is 
natural to consider a classical ground state as a refer- 
ence configuration. When one quantizes the spins along 
a local quantization axis in the direction of the spin ori- 
entation in the classical ground state, an interesting al- 
gorithm with open string-clusters emerges. The spins in 
each cluster are in the reference configuration and hence 




FIG. 1: kagome lattice (left) and frustrated square (or 
anisotropic triangular) lattice (right) consisting of three sub- 
lattices A, B, C . 

these clusters are protected from becoming unphysically 
large. However, the meron-concept does not apply to the 
open string-clusters, i.e. when these clusters are flipped, 
they are not independent but affect each other in their 
effect on the sign. Remarkably, one can still integrate 
out the spins analytically. This glues the open string- 
clusters together to the closed loop-clusters of the algo- 
rithm discussed before. While typical closed loop-clusters 
are hence larger than the correlation length correspond- 
ing to the classical order, they still represent physical 
correlated regions. In fact, they grow up to the length 
scale at which the signs, which are a manifestation of 
quantum entanglement, decorrelate. 

Even if the typical cluster-size is moderate, the inner 
Monte Carlo algorithm may not lead to an efficient can- 
cellation of signs. For example, there are cases in which 
the improved estimator (Sign}^ is not positive. Still, if 
such cases are rare, the sign problem is substantially re- 
duced. In order to optimize the performance of the al- 
gorithm, the numerical effort invested in the inner and 
outer Monte Carlo procedures must be properly balanced 
against each other. It pays off to invest a larger number 
of inner Monte Carlo sweeps on the larger sets Ac- In 
any case, the efficiency of the nested cluster algorithm 
must be investigated on a case by case basis. 

We now consider the Heisenberg antiferromagnet with 
uniform nearest-neighbor coupling J xy = J on the lat- 
tices illustrated in figure 1. The frustrated square lattice 
has an additional coupling J' along the diagonals. We 
have simulated large kagome lattices with up to V w 1000 
spins at moderate temperatures with (3 J » 1. Figure 2 
shows the probability distribution of the improved es- 
timator (Sign}^. Although sometimes it is negative, it 
still leads to an accurate determination of the average 
sign. We consider M s with z x = 1,-1,0 on sub-lattice 
A, B, C, respectively, which may signal coplanar spin or- 
der. As shown in figure 3, with increasing volume V both 
(Sign) + and (M 2 Sign} + decrease dramatically over nu- 
merous orders of magnitude, but are still accurately ac- 
counted for by the nested cluster algorithm. For example, 
with V = 882 spins (Sign)+ = 2.09(8) x 10~ 14 . A brute 
force approach would require an astronomical statistics of 



4 



1 1 1 1 


I 






pj = 1 




- 




V = 576 spins 




i iniiiiillllll 





-l.Oxlo'' -5.0X10" 10 00 5-OxlO" 10 I.OX10" 5 

< Sign > 



FIG. 2: Probability distribution of {Sign} i for the kagome lat- 
tice with V = 576 spins and ft J = 1. 
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FIG. 3: Volume- dependence of (Sign)+ and (MsSign) + 
(rescaled by 10~ 6 ) for the kagome lattice with (3 J = 1. 
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about 10 30 sweeps in order to achieve a similar precision. 
Figure 4 shows the coplanar staggered susceptibility Xs 
compared to the collinear Neel susceptibility xn- On the 
square lattice, frustration reduces the Neel order, while 
(at least for J' — J/4) the coplanar order is as weak as 
on the kagome lattice (and practically indistinguishable 
from it in figure 4). 

To conclude, in contrast to other Monte Carlo meth- 
ods, the nested cluster algorithm is capable of eliminating 
very severe sign problems for large systems, at least at 
moderate temperatures. This is useful, for example, for 
determining the couplings of frustrated magnets by com- 
parison with experimental finite temperature data. As 
we have demonstrated, although the nested cluster algo- 
rithm cannot reach very low temperatures, by studying 
appropriate susceptibilities one may still obtain valuable 
insights concerning possible types of order. Applications 
to frustrated antiferromagnets on various lattice geome- 
tries are currently in progress. 
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FIG. 4: Coplanar staggered susceptibility Xs and collinear 
Neel susceptibility xn as functions of the space-time volume 
(SV for the kagome as well as the frustrated (J' = J/4) and 
unfrustrated (J' = 0) square lattice at fixed space/time aspect 
ratio VV/PJ = 20. 
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